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Abstract Land surface phenology is widely retrieved from 
satellite observations at regional and global scales, and its 
long-term record has been demonstrated to be a valuable tool 
for reconstructing past climate variations, monitoring the dy- 
namics of terrestrial ecosystems in response to climate im- 
pacts, and predicting biological responses to future climate 
scenarios. This study detected global land surface phenology 
from the advanced very high resolution radiometer (AVHRR) 
and the Moderate Resolution Imaging Spectroradiometer 
(MODIS) data from 1982 to 2010. Based on daily enhanced 
vegetation index at a spatial resolution of 0.05 degrees, we 
simulated the seasonal vegetative trajectory for each individ- 
ual pixel using piecewise logistic models, which was then 
used to detect the onset of greenness increase (OGI) and the 
length of vegetation growing season (GSL). Further, both 
overall interannual variations and pixel-based trends were 
examined across Koeppen’s climate regions for the periods 
of 1982-1999 and 2000-2010, respectively. The results show 
that OGI and GSL varied considerably during 1982-2010 
across the globe. Generally, the interannual variation could 
be more than a month in precipitation-controlled tropical and 
dry climates while it was mainly less than 15 days in 
temperature-controlled temperate, cold, and polar climates. 
OGI, overall, shifted early, and GSL was prolonged from 
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1982 to 2010 in most climate regions in North America and 
Asia while the consistently significant trends only occurred in 
cold climate and polar climate in North America. The overall 
trends in Europe were generally insignificant. Over South 
America, late OGI was consistent (particularly from 1982 to 
1 999) while either positive or negative GSL trends in a climate 
region were mostly reversed between the periods of 1982— 
1999 and 2000-2010. In the Northern Hemisphere of Africa, 
OGI trends were mostly insignificant, but prolonged GSL was 
evident over individual climate regions during the last 3 
decades. OGI mainly showed late trends in the Southern 
Hemisphere of Africa while GSL was reversed from reduced 
GSL trends (1982-1999) to prolonged trends (2000-2010). In 
Australia, GSL exhibited considerable interannua! variation, 
but the consistent trend lacked presence in most regions. 
Finally, the proportion of pixels with significant trends was 
less than 1 % in most of climate regions although it could be as 
large as 10 %. 

Keywords Long-term global phenology • Interannual 
variation * Trend • Remote sensing 


Introduction 

Vegetation phenology reflects a complex suite of interactions 
among atmospheric, biospheric, and soil biogeochemical 
properties. Because vegetation phenology effectively reflects 
climate change, the Intergovernmental Panel on Climate 
Change (IPCC) reported that phenology is one of the simplest 
and most effective indicators of climate change (IPCC 2007). 
Long-term observations and recording of changes in plant 
phenology support efforts to understand trends in regional 
and global climate change, to reconstruct past climate varia- 
tions, to explore the magnitude of climate change impacts on 
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vegetation growth, and to predict biological responses to 
future climate scenarios. 

Vegetation phenology for a specific species is widely ob- 
served in fields. Crop calendars have been used in agriculture 
for thousands of years in China (Zhu and Yuan 1963), and 
phenological records have a history going back to the early 
1700 s in Europe (e.g., Sparks and Carey 1995) and to the 
1800 s in Japan (Lauscher 1978). Recently, several networks 
of field phenological observations have been established 
worldwide to record phenological events distributed spatially 
across a regional area. These networks include PI ant Watch in 
Canada, the National Phenology Network (NPN) in the USA 
(United States of America), the European Phenology 
Network, the Japan Phenological Eyes Network (PEN), and 
the UK (United Kingdom) Phenology Network. These net- 
works record the timing of species-specific events. More 
recently, a network of “phenocams” has been established in 
North America (NA), which currently includes more than 50 
sites (Richardson et al. 2009). The phenocams observe the 
seasonality of plant canopy in the field from digital webcams, 
which provide “near-surface” remotely sensed vegetation in- 
dices and have the potential to provide a rich source of 
quantitative data related to phenology at daily time scales. 

Throughout the last 3 decades, satellite remote sensing has 
become a widely used mechanism for monitoring phenology 
in vegetation communities, which is referred to as land surface 
phenoloy (LSP) (de Beurs and Henebry 2004; Friedl et al. 
2006). To detect long-term ecosystem variation and climate 
change at regional and global scales, LSP has been commonly 
retrieved from time series of the normalized difference vege- 
tation index (NDVI) from the advanced very high resolution 
radiometer (AVHRR), the VEGETATION instrument on- 
board the SPOT spacecraft, and the Moderate Resolution 
Imaging Spectroradiometer (MODIS) onboard NaSA’s 
Terra and Aqua spacecraft. AVHRR dataset, available since 
July 1981, provides the longest daily global observations and 
is commonly used for detecting land surface change at local 
and global scales. Different calibration methods generate sev- 
eral AVHRR NDVI datasets. NASA’s Global Inventory 
Modeling and Mapping Studies (GIMMS) group produced a 
global NDVI dataset at a resolution of biweekly 8 km from 
1981 to 2011 (Tucker et al. 2005). NASA/NOAA Pathfinder 
AVHRR Land (PAL) provided a global 10-day composite of 
8-km NDVI, which spans from July 1981 to December 1999 
(James and Kalluri 1 994). The USGS (US Geological Survey) 
EROS Data Center generated a 10-day composite of l -km 
AVHRR NDVI for the contiguous USA since 1989 
(Eidenshink 1992). NOAA NESDIS (National 
Environmental Satellite, Data, and Information Service) de- 
veloped a weekly 4-km global vegetation index (GVIx) prod- 
uct from 1981 to present (Kogan et al. 2003). The latest 
dataset is the Long Term Data Record (LTDR), which is 
designed to produce a consistent long-term record from the 


AVHRR and MODIS sensors (Pedelty et al. 2007). The 
AVHRR LTDR provides daily 0.05° (~5 km) reflectance data 
from 1981 to 1999. 

The time series of NDVI data has been extensively used for 
detecting phenological trend and interannual variation. 
AVHRR NDVI data have been commonly applied to retrieve 
the start and end of vegetation growth using a number of 
methods including NDVI threshold (Lloyd 1990; Fisher 
1994; Myneni et al. 1997), middle point (White et al. 1997), 
moving average (Reed et al. 1994), growing-degree-day- 
based NDVI model (de Beurs and Henebry 2004), Gaussian 
function (Jonsson and Eklundh 2002), and curvature change 
rate (Zhang et al. 2007). Analyzing long-term AVHRR- 
derived phenology has revealed an early trend in spring events 
and a late trend at the end of the growing season in middle- 
high latitudes (e.g., Myneni et al. 1997; Tucker et al. 2001; 
Zhou et al. 2001; de Beurs and Henebry 2005a; Chen et al. 
2005; Delbart et al. 2005; Piao et al. 2006; Zhang et al. 2007; 
Julien and Sobrino 2009; de Jong et al. 2011). The shift rate 
varies with different AVHRR time periods and research ap- 
proaches, ranging from 0.1 to 0.8 days/year in spring events 
and 0.1 to 1.1 days/year in growing season length. However, 
the late spring trend has also been found in southern NA from 
1982 to 2005 (Zhang et al. 2007) and in boreal forests from 
1993 to 2004 (Delbart et al. 2005). In contrast, it has been 
found that spring vegetation growth generally shows no sig- 
nificant trend from 1982 to 2006 over NA (Reed 2006; White 
et al. 2009). 

This study investigated interannual LSP variation using a 
daily enhanced vegetation index which was calculated from 
the AVHRR LTDR and the MODIS Climate Modeling Grid 
(CMG) records during 1981-2010. The temporal trajectory of 
the vegetation index was simulated using piecewise logistic 
models for land surface phenology detection after satellite 
data quality was explicitly employed. Further, phenological 
trend and interannual variation in both greenup onset and 
growing season length (GSL) were examined across 
Koeppen’s climate regions in each continent, respectively. 


Methodology 

Long-term daily vegetation index 

The enhanced vegetation index (EVI), calculated from reflec- 
tance in blue, red, and near-infrared bands, has advantages 
over NDVI in quantifying vegetation activity (Huete et al. 
2002). It reduces sensitivity to soil, nonphotosynthetically 
active vegetation (i.e., litter and woody tissues), and atmo- 
spheric effects but still remains sensitive to variation in canopy 
density where NDVI becomes saturated (Huete et al. 2002; 
Rocha et al. 2008). 


£) Springer 




Int J Biometeorol 


Recently, a two-band EVI (EVI2) has been proposed, 
which is calculated from red and near-infrared reflectances 
(Huete et al. 2006; Jiang et al. 2008). EV12 remains function- 
ally equivalent to EVI but can be derived from satellite sensors 
without blue band, such as AVHRR. Thus, EVI2 could be 
used :o monitor vegetation phenology and activity across a 
variety of ecosystems over a long-time period (Rocha and 
Shaver 2009). 

A long-term dataset of global EVI2 has been generated 
from the daily land surface reflectance in the AVHRR LTDR 
and the MODIS CMG records. The AVHRR LTDR (1981- 
1999' was produced by a NASA-funded REASoN project 
(Pedelty et al. 2007). It aims to produce a multidecade time 
series of reliable data from several sensors during the 1 980s 
and 1990s, compatible with data retrieved from MODIS. The 
LTDR (version 3) processes Global Area Coverage (G AC) data 
using the following approaches: (1) radiometric corrections 
were processed using ocean-clouds vicarious calibration algo- 
rithm (Vermote and Kaufman 1995), in which the independent- 
ly derived sets of AVHRR calibration coefficients were consis- 
tent with an accuracy of 1 %; (2) viewing and illumination were 
adjusted using Bidirectional Reflectance Distribution Function 
(BRDF) techniques; (3) rigorous cloud and cloud shadow were 
screened using Cloud Advanced Very High Resolution 
Radiometer (CLAVR-1) (Stowe et al. 1999); and (4) water- 
vapor observation was corrected using reanalysis data from 
the NOAA NCER As a result, the LTDR dataset (version 3) 
consists of daily data including information from AVHRR 
channels 1 to 5, as well as ancillary data describing sun- 
sensor-target geometry at a spatial resolution of 0.05° (~5 km), 
which is the same as MODIS Climate Modeling Grid (CMG) 
products (since 2000) (http://ltdr.nascom.nasa.gov/cgi-bin/ltdr/ 
ltdrPage.cgi?fileNameHMODlS2LTDR). 

A long-term data record of spectral vegetation index (VI, 
which are NDVI and EVI2) from both AVHRR LTDR and 
MODIS CMG was later generated with the support of a 
NASA MEaSUREs project. To merge VI from multiple- 
sensor data, data continuity and compatibility have been sci- 
entifically investigated. Analysis of AVHRR GAC broader 
band observations and Terra MODIS CMG narrower bands, 
which are simulated from Hyperion data, has demonstrated 
capability in combining EVI2/NDVI from different sensors 
(Kim et al. 2010). Because NOAA-14 AVHRR from 2000 to 
2001 is not recommended for scientific uses due to significant 
orbital drift, there is no overlapping period of observations 
between AVHRR LTDR and MODIS CMG. Thus, SPOT-4 
VGT was applied to bridge these two datasets (Tsend-Ayush 
et al. 2012). The long-term daily VI data from AVHRR and 
MODIS data at a spatial resolution of 0.05° were generated 
using a set of functions with two linear polynomials (Tsend- 
Ayush et al. 2012; Yoshioka et al. 2012). The VI relationship 
equations eliminate the soil brightness effect by utilizing the 
vegetation isoline equations and reduce the effects from 


biophysical properties by separating land cover types 
(Yoshioka et al. 2012). We obtained the long-term daily 
EVI2 over the last 30 years from the University of Arizona 
(http://vip.arizona.edu/viplab_data explorer). 

Detection of vegetation phenology 

To establish a long-term record of global LSP from satellite 
data, phenological metrics are defined according to vegetation 
seasonal growing cycles. Briefly, a seasonal cycle of vegeta- 
tion growth consists of a greenup phase, a maturity phase, a 
senescent phase, and a dormant phase (Zhang et al. 2003). 
These four phases are characterized using four phenological 
transition dates in the time series of greenness data: (1) 
greenup onset, the date of onset of greenness increase 
(OGI); (2) maturity onset, the date of onset of greenness 
maximum; (3) senescence onset, the date of onset of green- 
ness decrease; and (4) dormancy onset, the date of onset of 
greenness minimum. Furthermore, GSL is defined as the time 
difference between greenup onset and dormancy onset. 

To identify the global phenological metrics, we used a 
piecewise-logistic-model-based LSP detection algorithm 
(PLM-LSPD, Zhang et al. 2003). This PLM-LSPD has been 
demonstrated to be effective in depicting the seasonality of 
satellite greenness across various ecosystems (Zhang et al. 
2003, 2006). This was then applied to investigate seasonal 
vegetation growth using webcam data (Richardson et al. 2006; 
Kovalskyy et al. 2012), Landsat TM data (e.g., Fisher et al. 
2006; Kovalskyy et al. 2012), AVHRR data (e.g., Zhang et al. 
2007; Julien and Sobrino 2009), and MODIS data (e.g., Zhang 
et al. 2003, 2006; Ahl et al. 2006; Liang et al. 2011). Thus, the 
PLM-LSPD is believed to be applicable for the detection of 
long-term global vegetation phenology from EVI2 data. 

Here, we briefly introduce the PLM-LSPD for processing 
temporal EVI2 data for phenology detection. First, we 
smoothed the time series of EVI2 based on the following 
assumptions: (1) vegetation growth is a continuous process 
without sharp increases or decreases, (2) the time between 
successive peaks in separate vegetation growth cycles is lon- 
ger than 3 months in forests and longer than 2 months for other 
plant functional types, (3) the magnitude of EV12 values is 
generally lowered by contamination from clouds and atmo- 
sphere, and (4) a single isolated maximum in EVI2 value in a 
time series is not always reliable. To reduce the data size and 
computation time, we first generated a 3-day composite 
dataset from daily EV12 data by selecting the best quality data 
and using the maximum value composite. We then removed 
the observations of clouds and snow cover from temporal 
EVI2 data. The cloud and snow covers are described in the 
quality assurance flag in the LTDR and MODIS CMG 
datasets. Additionally, the background EVI2 value at each 
pixel, which represents the minimum EVI2 consisting of soil 
and seasonally stable vegetation in an annual time series 
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(Zhang et al. 2007), was identified and used to replace EV12 
values that varied irregularly during a winter period. The data 
gaps caused by clouds were filled by linear interpolation using 
neighboring good quality data. The time series of EVI2 data at 
each pixel was finally smoothed using a Savitzky-Golay filter 
and a running local median filter. 

EVI2 time series in a growing phase (EVI2 consistent 
increase) and a senescent phase (EVI2 consistent decrease) 
was modeled using a logistic function, respectively. 
Phenological transition dates within each growth or senes- 
cence phase were identified using the rate of change in the 
curvature of the modeled curves. Specifically, transition dates 
correspond to the day of year (DOY) on which the rate of 
change in curvature in the EVI2 data exhibits local minima or 
maxima (Zhang et al. 2003). Note that if the annual variation 
in EVI2 was subtle (<0.04), such as in some evergreen forests 
and shrublands, phenological metrics were not retrieved. 
Moreover, if the annual EVI2 time series was very flat in 
some shrublands, the curvature change rate might produce a 
very early timing of greenup onset. In this case, a 5 % thresh- 
old was used to determine the start and the end of a growing 
season. 

GSL was calculated using the difference between the 
timing of greenup onset and dormancy onset. This was 
straightforward in the regions where vegetation generally 
grew from spring to autumn within a calendar year. For the 
regions where a growing cycle spanned 2 different calendar 
years, GSL was assigned to the year that greenup onset 
occurred. 

Determination of the quality of an annual time series EVI2 

The accuracy of phenology detection is strongly dependent on 
the quality of time series EV12 data. A smooth approach could 
reduce some of the uncertainties, but missing observations in 
the time series greatly reduce the detection precision of 
phenologic metrics (Zhang et al. 2009). We used the propor- 
tion of good quality observations ( Pqa ) during a vegetation 
growing season (the time period between greenup onset and 
dormancy onset) to evaluate the phenology detections: 

Pqa = Nqaj T 

where T is the total number of 3-day EVI2 during a growing 
season, and Nqa is the number of good quality observations 
within a vegetation growing season. An EVI2 is counted as 
one good observation if there is one good value within a 
moving window of three 3-day EVI2. This is due to the fact 
that the error in phenology detection is small if the annual 
greenness cycle is simulated using a logistic model from EVI2 
data with the temporal resolution finer than 8 days (Zhang 
et al. 2009). 


Investigation of trend and interannual variation in phenology 

The trend and interannual variation in both OGI and GSL 
were analyzed across the globe. This analysis was performed 
based on climate regions respectively in each continent. This 
was due to the fact that the magnitude and rate of climate 
change primarily vary with climate regions which are com- 
monly classified based on monthly temperature and precipi- 
tation and that climate change can be effectively indicated by 
the variation in vegetation phenology. In other words, if the 
climate variables change within a climate region, vegetation 
phenology would vary correspondingly, although phenologi- 
cal timing could vary with vegetation types. To examine 
phenological variation in a regional or local area, ecoregion 
or phenoregion could provide localized patterns because 
ecoregions distinguish areas that share common climatic char- 
acteristics and then vegetation properties (Olson et al. 2001), 
and phenoregions characterize seasonal variation in NDVI 
and climatology (White et al. 2005; Hargrove et al. 2009). 
However, in investigating the overall variation in long-term 
phenology across the globe, it is assumed that the phenolog- 
ical variation in climate regions across various geographic 
areas could reflect the general pattern of climate change. 
Therefore, we first stratified both OGI and GSL by geographic 
regions and climatic classes, respectively (Fig. 1). The geo- 
graphic regions were divided into North America (NA), South 
America (SA), Europe (EU), Asia (AS), Australia (AU), 
Northern Hemisphere of Africa (NAF), and Southern 
Hemisphere of Africa (SAF). The climatic classes, selected 
from Koeppen’s climate classification, were developed by the 
FAO-SDRN Agrometeorology group (http://www.fao.org/nr/ 
climpag/globgrids/KC classification en.asp, 2006) (Fig. 1). 
This system consists of the following classes: (1) tropical 
climate — fully humid season (Af), monsoon type (Am), dry 
summer (As), and dry winter (AW); (2) dry climate— steppe 
climate (BS) and desert climate (BW); (3) temperate 
climate — fully humid season (Cf), dry winter (Cw), and dry 
summer (Cs); (4) cold climate — fully humid season (Df), dry 
winter (Dw), and dry summer (Ds); and (5) polar 
climate — tundra (ET) and frost climate (EF). 

The time series of OGI and GSL was finally applied to 
detect the trends for each climate region using a linear regres- 
sion model. The trends were calculated for the two periods of 
1982-1999 and 2000-2010, separately. This was due to the 
fact that phenology from 1982 to 1999 was derived from 
AVHRR data which were widely used for phenological trend 
detection in previous studies. In contrast, phenology from 
2000 to 2010 was detected from MOD1S data, which have 
better quality because MODIS land spectral bands were ex- 
plicitly designed for land surface monitoring (Justice et al. 
1997). Moreover, we found that the inconsistency between 
AVHRR LTDR and MODIS CMG still remains in some 
regions although great efforts have been exerted to continue 
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Fig. 1 Koeppeirs climate classification (based on FAO-SDRN Agrometeorology group, 2006) and geographic regions 


these two datasets. Such sensor inconsistency could signifi- 
cantly affect trend analysis (de Beurs and Henebry 2005b). 

In each time period, we calculated overall trend and pixel- 
based trend in a climate region, separately, after the geograph- 
ic projection (0.05°) was reprojected to the Goode homolosine 
projection (an equal-area project) at a spatial resolution of 
5 km. The overall trend was calculated from a time series of 
average phenological values in a region. The pixel-based trend 
was used to determine the pixels with significantly positive 
and negative trends and the corresponding proportion of 
pixels (area) in a climate region. Crop land was excluded in 
this analysis because crop phenology is controlled by human 
activities. This data processing allowed us to understand well 
both the overall changes and the percentage of changing 
pixels that were significantly affected by climate change with- 
in a climate region. 


Results 

Quality of annual time series of EVI2 

Reliable phenology retrieval is strongly dependent on the 
quality of annual time series of satellite observations. 
Figure 2 presents the spatial pattern of the average percentage 
of good quality observations. Overall, the percentage of good 
observations was much higher in MODIS observations 
(2000-2010) than in AVHRR obseivations (1982-1999). In 
MODIS observations, the percentage of good observations 
during a vegetation growing season was generally larger than 
75 % except in tropical climates, England, and northwestern 
NA. In contrast, good observations from AVHRR data were 
less than 75 % in most regions except for dry climates. 
However, good observations were less than 50 % in large 
parts of Western Europe, East Asia, the Amazon basin, and 
West Africa (Fig. 2). This is due to the fact that MODIS CMG 
data were aggregated from 500-m pixels with high quality. In 
contrast, AVHRR LTDR was generated from 4-km pixels 


which were sampled onboard the satellite by selecting every 
third 1-km local area coverage (LAC) data line, four of every 
five samples along this line, and averaging these data. 

Interannual variation in the onset of greenness increase 

Figure 3 shows the global climatology of the onset of green- 
ness increase from 1982 to 2010. The spatial pattern in OGI 
reflects both broad-scale variations in controlling mechanisms 
related to climate and more local factors related to land cover 
and human activities (Fig. 3a). Specifically, OGI spatially 
shifted along latitude from 30°N northward. Zonal patterns 
indicate that OGI varied at a rate of 2-3 days per degree of 
latitude in NA, Europe, and Asia. This rate of northward 
progression is similar to Hopkins’ law (1938; Reader et al. 
1974). In NA, OGI shifted from late February in southern 
USA, April in northern USA, and to the end of June in the 
northern tundra. The dependence on latitude was spatially 
interrupted by elevation and human activities. The elevation 
effect was pronounced in large mountains such as the 
Himalayas in Asia, the Alps in Europe, and the Appalachians 
in America. Moreover, OGI in agricultural lands was frequent- 
ly distinct from that of surrounding natural vegetation because 
of the controls applied by human management. This was 
highly evident in central NA, where OGI occurred much later 
in the Mississippi River valley and the midwestem agricultural 
heartland, relative to the surrounding natural vegetation. 

OGI pattern was generally complex and irregular in dry 
climate (arid and semiarid regions) because of the highly 
spatial variability in water availability. However, regular pat- 
terns in local regions were commonly found. For example, in 
the Northern Hemisphere of Africa (the Sahelian and sub- 
Sahelian regions), OGI shifted smoothly from early March at 
around 6.5°N to mid-October in the boundary between the 
Sahel and the Sahara desert (17.9°N). The shift rate was about 
20 days per degree of latitude; this reflected the start of the 
rainy season which was controlled by the migration of the 
Intertropical Convergence Zone (ITCZ) (Zhang et al. 2005). 
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Fig. 2 Average percentages of 
good observations during a 
vegetation growing season in 
1980s (a), 1990s (b), and 2000s 
(c). Note that no percentage was 
calculated for the pixels if no 
phenoiogy was detected for any 
years during each period 



b 





In southern South America, green leaves emerged in the 
boreal summer and gradually pushed northward at a rate of 
about 3 days per latitude from 45° to 20°S. Along the 
Brazilian Highlands (in the direction from 60°W and 39°S to 
35°W and 5°S), OGI shifted from July to the following 
February at a rate of about 0. 12 days/km. 

Figure 3b shows the spatial pattern in standard deviation 
(STD) of temporal OGI from 1982 to 2010. It reveals that 
STD was less than 1 5 days in most regions of middle-high 
latitude across the Northern Hemisphere, particularly across 
Eurasia. In contrast, STD could be larger than 40 days in 
tropical and subtropical regions. This is likely associated with 
climate forces in controlling vegetation phenology. The vari- 
ance of temperature-controlled phenological timing was rela- 
tively small because seasonal temperature is generally regular 
among different years. However, rainfall seasonality could 
have large variances among different years (Lotsch et al. 
2003; Zhang et al. 2010), which results in large OGI STD in 
water-controlled arid climate regions. In tropical rainforests, 
such as the Amazon, vegetation greenness phenology is 


affected by sunlight during the dry season (Huete et al. 
2006), which leads to a large interannual variation in OGI. 

Figure 4 presents interannual variation in OGI across dif- 
ferent geographic and climate regions. Generally, OGI pre- 
sented considerable intcrannual variability in tropical, temper- 
ate, and diy climates while it was relatively consistent in cold 
and polar climates. The STD of interannual OGI was gener- 
ally less than 8 days in cold and polar climates and 1 0 days in 
temperate climate over the Northern Hemisphere. It was less 
than 21 days in tropical seasonally dry climates over the 
Southern Hemisphere. In diy climate, the STD was as large 
as 1 4 days in western N A, 7 days in Eurasia, 23 days in Africa, 
and 1 7 days in South America. 

The overall trends in OGI showed complexity during the 
past 3 decades (Fig. 4 and Tables 1 and 2). In the northern 
tundra (ET), OGI, overall, showed an early trend at a rate of 
0.25 days/year from 1982 to 1999 and 0.68 days/year from 
2000 to 2010 in NA, but the early trend was not significant in 
both Europe and Asia. Accordingly, trend analysis of individ- 
ual pixels in ET showed that the proportion of pixels was 2.2- 
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Fig. 3 Spatial pattern in OGI 
climatology based on data from 
1982 to 2010. a Average of OGI 
and b standard deviation. Note 
that the white color indicates no 
phenology detection because of 
either permanent snow cover or 
desert 


a 





b 




3.6 °/c for significantly early trends and less than 0.2 % for late 
trend during these two periods in NA. It was less than 0.5 % 
for late trend and 1 .2-1 .5 % for early trend in Asia, while the 
proportion was less than 1 % in either directions over Europe 
(Figs. 5 and 6). 

OGI tended to be early in the cold climate of Asia and NA 
(Table 1). In Df climate, OGI over NA showed a significantly 
early shift (0.38 days/year) from 1982 to 1999 but insignifi- 
cantly early shift from 2000 to 2010. No significant trend was 
found in either Europe or Asia. In Dw climate, which covers a 
large area of Asia, OGI came earlier at a significant rate of 
0.22 days/year from 1982 to 1999. In contrast, OGI presented 
consistent early trends in Ds climate across the Northern 
Hemisphere although the trend was insignificant in Europe. 
Pixel-based trends revealed that the proportion of pixels with 
significantly early shifts was 10 and 6 % of the Df climate in 
NA from 1 982 to 1 999 and from 2000 to 2010, respectively, 4 
and 2 % in Europe, and 5 and 3 % in Asia (Figs. 5 and 6). 
Accordingly, the area with significantly late shift was 1 and 
2 % in NA, 5 and 7 % in Europe, and 1 and 3 % in Asia. In Dw 
climate, the proportion was 4 % for early shift and less than 
1 % for late shift from 1982 to 1 999, and it was 1-2 % in both 
directions from 2000 to 2010 in Asia. The proportion for 
significant trends was less than 0.1 % in NA and Europe, 
respectively. 

In temperate climate, OGI showed slightly early trends in 
Cs from 1 982 to 2010 while it became late in Cw from 2000 to 
2010 over Asia. However, it presented no significant trends in 
both NA and Europe. At pixel level, area with either early or 


late trend was basically less than 1 % across the Northern 
Hemisphere. 

In dry climate, the significant trend in OGI only occurred in 
BW of Asia (Table 1). It became early at a rate of 0.83 days/ 
year from 1982 to 1999. The proportion of pixels with signif- 
icant shift was less than 1 % in various regions except for BS 
climate (Figs. 5 and 6). In the latter, the proportion was 1.1 % 
(late shift from 2000 to 2010) in NA, 1.1 % (late shift from 
1982 to 2010), and 1.3 % (early shift from 2000 to 2010) in 
Europe. 

In South America, OGI mostly presented late trends during 
the last 3 decades (Fig. 4 and Table 2). The significantly late 
trend occurred in all climates from 1 982 to 1 999 except for Cs 
and ET while it only occurred in Cw from 2000 to 2010. The 
proportion of pixels with significantly late shift could be as 
large as 5 % in an individual climate region, which was much 
larger than that for early shift (Figs. 5 and 6). 

Overall, OGI showed no significant trend in the Northern 
Hemisphere of Africa although early trends were transited to 
late trends from 1982-1999 to 2000-2010 (Table 2). 
However, the proportion of pixels could be as large as 3 % 
in a climate region, and it was generally larger for early shift 
than that for late shift (Figs. 5 and 6). Following the climate 
types, which are BW, BS, AW, and As, along latitude from 
northern to southern regions, the proportion for significantly 
early OGI shifts from 2000 to 2010 increased from 0.7 to 
3.2 %, gradually. In contrast, OGI in the Southern Hemisphere 
of Africa mainly showed late trends. The significant trend 
appeared in BS, Bw, Cw from 1982 to 1999, and Cf from 
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Fig. 4 Interannual variation in mean OGI across various geographic and climate regions 

2000 to 2010 (Table 2). The proportion of pixels vtas much based trend from 2000 to 20 10 was about 4 days/year in either 

larger for late shift than that for early shift (Figs. 5 and 6). early or late shift with an area of 1-4.5 % in AW, BS, BW, and 

In Australia, significantly early OGI trend only appeared in Cf, respectively (Figs. 5 and 6). Similar pattern appeared 

tropical climate ( AW, 3. 2 days/year) from 2000 to 20 10. Pixel- during the period from 1982 to 1999. Noticeably, the 


Table 1 Overall trend in the onset of greenness increase (days/year) in North America (NA), Europe (EU), and Asia (AS) 


OGI 


AW 

BS 

BW 

Cf 

Cw 

Cs 

Df 

Dw 

Ds 

ET 

NA 

1982-1999 

N/A 

0.22 

-0.56 

0.01 

-0.36 

0.08 

-0.38 

0.02 

-0.45 

-0.25 


2000-2010 

N/A 

1.20 

-1.73 

0.11 

-1.03 

0.62 

-0.41 

0.05 

-1.21 

-0.68 

EU 

1982-1999 

N/A 

-0.02 

-0.35 

-0.2 

N/A 

0.21 

0.02 

N/A 

-0.17 

0.13 


2000-2010 

N/A 

-0.13 

-0.28 

-0.31 

N/A 

-0.68 

0.02 

N/A 

-0.61 

-0.6 

AS 

1982-1999 

-0.32 

0.01 

-0.83 

-0.28 

-0.11 

-0.41 

-0.21 

-0.22 

-0.27 

-0.15 


2000-2010 

1.24 

-0.38 

0.41 

-0.3 

0.87 

-0.95 

0.18 

0.06 

-0.39 

-02 


The bold value indicates significant trend ( P value<0.1) 

N/A net applicable, OGI onset of greenness increase, AW tropical climate — diy winter, BS dry climate — steppe climate, BW dry climate — desert climate, 
Cf temperate climate — fully humid season, Cw temperate climate — dry winter, Cs temperate climate — diy summer, Df cold climate — fully humid 
season, Dw cold climate-dry winter, Ds cold climate — dry summer, ET polar climate — tundra 
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Table 2 Overall trend in the onset of greenness increase (days/year) in South America (SA), the Northern Hemisphere of Africa (NAF), the Southern 
Hemisphere of Africa (SAF), and Australia (AU) 


OGI 


As 

AW 

BS 

BW 

Cf 

Cw 

Cs 

ET 

SA 

1982-1999 

3.13 

2.24 

1.51 

0.71 

1.14 

1.42 

-0.26 

0.24 


2000-2010 

0.72 

1.09 

0.82 

0.57 

-0.38 

2.21 

1 

0.33 

NAF 

1982-1999 

0.07 

-0.36 

-0.24 

-0.2 

-0.44 

-0.11 

-0.77 

N/A 


2000-2010 

-0.67 

0.05 

0,19 

0.55 

0.63 

0.93 

-0.65 

N/A 

SAF 

1982-1999 

0.02 

1.16 

1.63 

0.13 

0.29 

138 

-0.18 

N/A 


2000-2010 

1.22 

0.55 

0.78 

1.37 

234 

-0.09 

0.07 

N/A 

AU 

1982-1999 

N/A 

0.97 

1.19 

-0.88 

0.41 

-3.00 

-0.10 

N/A 


2000-2010 

N/A 

-3.18 

-2.11 

-0.73 

0.71 

-3.16 

-3.45 

N/A 


The bold value indicates significant trend (. P value<0.1) 

N/A not applicable, OGI onset of greenness increase. As tropical climate — dry summer, AW tropical climate — dry winter, BS dry climate — steppe 
climate, BW dry climate — desert climate, O' temperate climate — fully humid season, Cw temperate climate — dry winter, Cs temperate climate — dry 
summer, ET polar climate — tundra 
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proportion for early shift was consistently large (~4 %) in BS 
from 1982 to 2010. Overall, OG1 tended to appear earlier 
across the region. 

Spatial pattern and interannual variation in growing season 
length 

Figure 7 shows the spatial pattern of average GSL in the 
1980s, 1990s, and 2000s, respectively. As expected, the pat- 
tern is very similar among the three periods. In the Northern 
Hemisphere, GSL showed a clear regular gradient along the 
latitude. It varied from about 50 days in the tundra area to 
around 300 days in temperate climate regions. Similar to the 
variation in OGI, the regular pattern in GSL was interrupted 
by human activity and elevation. Indeed, the spatial shift of 
GSL is also dependent on elevation. GSL varied from 
220 days in Switzerland (north of the Alps) to about 90 days 


in the Alps Mountains and to 240 days in Italy (south of the 
Alps). 

In the Southern Hemisphere, the spatial pattern in GSL was 
complex. Regular patterns appeared in local regions. For 
example, GSL varied from 125 days in Namibia to 260 days 
in Zimbabwe along a direction of southwest to northeast. It is 
worth noting that GSL in tropical forests could be shorter than 
that of other vegetation types because GSL in evergreen 
forests represents a cycle period from greenness increase to 
decrease. GSL was remarkably short in the horn of Africa 
because it only represented the first growing cycle in Fig. 7, 
where two vegetation growing cycles followed two rainy 
seasons (Zhang et al. 2005). 

Interannual variation in GSL showed a considerable spatial 
difference across various climate regions. It was relatively 
small in cold and polar climate regions with a STD of gener- 
ally less than 15 days. In contrast, the STD was larger than 
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Fig. 7 Spatial pattern in average GSL (a, b, and c) and interannual standard deviation (d, e, and f). Note the first growing season is presented if there are 
multiple seasons within a year. Gray indicates no phenology detection 


20 days in dry and tropical climate regimes. It was largest in 
tropical forests, where it could be more than 40 days. 

Long-term trend in GSL 

The GSL trend varies considerably with geographic and cli- 
mate regions (Figs. 8, 9, and 10, Tables 3 and 4). In NA, 
significantly prolonged GSL was most prominently observed 
from 1982 to 2010 although GSL became shorter insignifi- 
cantly in BS and Cs from 2000 to 2010 (Fig. 8 and Table 3). 
Particularly, GSL significantly increased consistently in Ds 
and ET during the past 30 years. This pattern was also 
reflected in pixel-based trends (Figs. 9 and 10). The propor- 
tion cf pixels was 1-3 % for prolonged GSL in ET and Ds 
while it was less than 0.4 % for shortened GSL. In Df, 
prolonged GSL occurred in 11 % of the area from 1982 to 
1999 and 4 % from 2000 to 2010, while shortened GSL 
appeared in 2 and 7 %, respectively. In other climate regions, 
the area that changed significantly was less than 1 %. 

GSL tended to be short in most climate regions across 
Europe from 1982 to 2010. It was significantly shortened in 
Cs from 1982 to 2010 and Ds from 2000 to 2010 (Table 3). 
Although the proportion of pixels was small (<1 %), there 
were more pixels with shortened GSL than those with 


prolonged GSL from 1982 to 2010 in ET, Ds, Cs, and BS, 
respectively. In Df, shortened GSL appeared in 4 % (1982— 
1999) and 10 % (2000-2010) of area while prolonged GSL 
occurred in 6 and 3 %, respectively (Figs. 9 and 10). 

Significantly prolonged GSL in Asia presented in BS 
(2000-2010), BW (2000-2010), and CS (1982-2010), while 
shortened GSL appeared in BS (1982-1999) and Df (2000- 
2010) (Table 3). The proportion of pixels with significant GSL 
shift was relatively large in Df and Dw while it was smaller 
than 1 % in other climate regions. In Df, it was 2 % ( 1 982— 
1999) and 3 % (2000-2010) for prolonged GSL, and 3 % 
(1982-1999) and 7 % (2000-2010) for shortened GSL, re- 
spectively. In Dw, the proportion was 3, 3, 1, and 1 %, 
respectively. 

In South America, GSL mainly showed opposite overall 
trends between 1982-1999 and 2000-2010 although the 
trends were not necessarily significant (Fig. 8 and Table 4). 
Significant opposite trends only occurred in Cf, where GSL 
increased (0.9 days/year) from 1982 to 1999 while it de- 
creased (1.7 days/year) from 2000 to 2010. Moreover, the 
GSL could transit from positive to negative trends across 
different climate regions. Indeed, GSL presented a prolonged 
trend in dry climate (BS and B W) while it showed a shortened 
trend in temperate climate (Cf and Cw) from 2000 to 2010. In 
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Fig. 8 Intemnnual variation in mean GSL across various geographic and climate regions 


general, there were more pixels with reduced GSL than 
increased GSL in AW, Cf, and Cw from 1982 to 2010 
(Figs. 9 and 10). 

GSL was generally prolonged across the Sahelian and sub- 
Sahelian regions where AW, BS, and BW are dominated, 
although the significant overall trend only occurred in BS 
from 1982 to 1999 (Fig. 8 and Table 4). However, GSL 
reduced from 1982 to 2010 in As with a significantly short- 
ened trend from 1982 to 1999. Overall, the proportion of 
pixels with prolonged GSL was larger than that with reduced 
GSL. It was 2-5 % for increased GSL in AW and BS from 
1982 to 2010, while the proportion was less than 1 % for 
decreased GSL (Figs. 9 and 10). In As and BW, there were 
slightly more pixels with increased GSL than decreased GSL 
from 2000 to 2010 while the proportion for either shift direc- 
tions was equivalent from 1982 to 1999. 

In the Southern Hemisphere of Africa, the significant GSL 
trend only appeared in As from 1982 to 1999. However, the 


proportion of pixels with significant trends could be as large as 
4 % in a climate region. The proportion for significantly 
reduced GSL was larger than that for prolonged GSL in AW, 
BS, BW, and Cw from 1982 to 1999, while this pattern was 
reversed from 2000 to 2010. As a result, the trend from 1982 
to 2010 was generally insignificant. 

In Australia, GSL presented an increased trend in BW 
(2.9 days/year) and a decreased trend in Cw (1.5 days/year) 
from 1982 to 1999. There were no significant trends across 
other climate regions during the two periods of 1982-1999 
and 2000-2010. However, the proportion of pixels showing 
negative trends was larger than that for positive trends in AW 
and BW from 2000 to 2010, which was similar in both 
directions from 1982 to 1999 (Figs. 9 and 10). In BS, the 
proportion for prolonged GSL was larger than that for reduced 
GSL from 1982 to 2010 although it was larger than 2 % in 
both positive and negative trends. In other climate regions, the 
area with significant trends was very small (<1 %). 
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Fig. 9 The shift rate averaged from significant pixel-based trends (P<0. 1 ) and the proportion of pixels with significant trends in growing seasonal length 
from 1982 to 1999 


Discussion 

A long-term global LSP dataset was produced from multiple- 
sensor satellite data in this study, however, to be aware of the 
limitation and quality that could improve our interpretation of 
the impacts of climate change on ecosystem dynamics. 
Although the abiotic noises caused by cloud cover and snow 
appearance were explicitly reduced in simulating temporal 
vegetative EV12 trajectory, the uncertainties in the reconstruct- 
ed EVI2 series are large if cloud cover is consistent during a 
vegetation growing season. Further, the annual EVI2 trajecto- 
ries may not well reflect subtle seasonal cycles in some 
evergreen forests and shrublands. Thus, the retrieved pheno- 
logical metrics are likely not so reliable as those from decid- 
uous forests, croplands, and grasslands. Moreover, temporal 
MOD1S data have a relatively higher quality than AVHRR 
data, and their continuity remains uncertain. As a result, these 


factors likely affect the phenological trend analysis from 1982 
to 2010. 

Validation is required to understand well the accuracy of 
long-term global LSP detections. It ensures the reliability of 
long-term trend and interannual variation of phenology in 
response to climate change. To validate the accuracy of phe- 
nology detection, sufficient field measurements comparable to 
a satellite footprint are needed. This requires the field data to 
reconcile with satellite-based phenological observations, 
which is currently extremely challenging. The validation ef- 
fort will become more practical with the inclusion of obser- 
vations from webcam (Richardson et al. 2009) and landscape 
measurements upscaled from field observations (Liang et al. 
2011). However, such data are currently very limited for 
validating our global LSP product. 

Comparison of the satellite-derived phenological trends 
from various studies could improve our understanding of the 
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Fig. 10 The shift rate averaged from significant pixel-based trends (P -0.1) and the proportion of pixels with significant trends in growing seasonal 
length from 2000 to 2010 


complexity in trend detections. The phenological trend and process the time series of satellite data. First, time series of 

interannual variation vary greatly with the methods used to NDVI and EVI2 are not temporally consistent, and EVI2 has 


Table 3 Overall trend in the growing season length (days/year) in North America (NA), Europe (EU), and Asia (AS) 


GSL 


AW 

BS 

BW 

Cf 

Cw 

Cs 

Df 

Dw 

Ds 

ET 

NA 

1982-1999 

N/A 

0.04 

1.14 

0.19 

0.87 

0.35 

0.38 

026 

0.43 

0.25 


2000-2010 

N/A 

-1.08 

0.38 

0.3 

1.27 

-0.37 

0.25 

0.65 

1.27 

0.65 

EU 

1982-1999 

N/A 

0.39 

-0.41 

-0.09 

N/A 

- 0.74 

0.14 

N/A 

0.16 

-0.02 


2000-2010 

N/A 

-1.33 

0.75 

0.46 

N/A 

- 1.86 

-0.31 

N/A 

-2.19 

0.55 

AS 

1982-1999 

0.02 

-0.46 

0.28 

-0.27 

0.39 

1.04 

0.05 

0.15 

0.12 

0.03 


2000-2010 

-0.86 

1.63 

0.57 

0.34 

-0.03 

3.75 

-0.44 

0.02 

0.25 

0.22 


The bold value indicatesignificant trends ( P value<0.1) 

N/A net applicable, GSL length of vegetation growing season, AW tropical climate — dry winter, BSdry climate — steppe climate, BWdry climate — desert 
climate, Cf temperate climate — fully humid season, Cw temperate climate — dry winter, Cs temperate climate — dry summer, Df cold climate — fully 
humid season, Dw cold climate — dry winter, Ds cold climate — dry summer, E T polar climate— tundra 
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Ikble 4 Overall trend in the growing season length (days/year) in South America (SA), the Northern Hemisphere of Africa (NAF), the Southern 
Hemisphere of Africa (SAF), and Australia (AU) 


GSL 


As 

AW 

BS 

BW 

Cf 

Cw 

Cs 

ET 

SA 

1982-1999 

0.89 

0.8 

-0.36 

0.15 

0.85 

0.01 

0.35 

0.3 


2000-2010 

-0.75 

-0.71 

0.87 

1.61 

- 1.67 

-239 

-0.41 

-0.52 

NAF 

1982-1999 

0.94 

0.18 

0.75 

0.66 

- 1.26 

0.59 

0.03 

N/A 


2000-2010 

-0.37 

0.36 

0.32 

-0.07 

-0.74 

- 1.62 

-0.45 

N/A 

SAF 

1982-1999 

0.96 

0.25 

-0.35 

0.26 

0.13 

0.03 

-0.26 

N/A 


2000-2010 

0.04 

0.23 

0.83 

-1.56 

-0.79 

0.06 

-0.87 

N/A 

AU 

1982-1999 

N/A 

-0.23 

0 

2.86 

0.28 

-134 

-0.15 

N/A 


2000-2010 

N/A 

-1.65 

-1.54 

-2.27 

0.08 

-0.18 

-0.85 

N/A 


The bold value indicates significant trend (, P value<0.1) 

N/A not applicable, GSL length of vegetation growing season. As tropical climate — dry summer, AW tropical climate— dry winter, BS dry climate — 
steppe climate, BW dry climate — desert climate, Cf temperate climate — fully humid season, Cw temperate climate — dry winter, Cs temperate climate — 
dry summer, ET polar climate — tundra 


advantages over NDVI in monitoring vegetation seasonal 
dynamics (Huete et al. 2006; Jiang et al. 2008; Rocha 
et al 2008). Second, although the same phenological term is 
used for defining a phenological event, such as OGI, the 
corresponding biophysical meaning is different among vari- 
ous methods such as absolute threshold (e.g., Lloyd 1990; 
Fisher 1994; Myneni et al. 1997; Zhou et al. 2001), relative 
threshold (e.g., White et al. 1997; Jonsson and Eklundh 2002; 
Delbart et al. 2005; Karlsen et al. 2006), quadratic model (de 
Beurs and Henebry 2005a), moving average (Reed et al. 
1994), change rate in NDVI (Moulin et al 1997; Tateishi 
and Ebata 2004; Piao et al. 2006), and curvature change rate 
(Zhang et al. 2003). Third, temporal greenness trajectories 
reflected seasonal vegetation dynamics are always poorly 
reconstructed although various smoothing approaches have 
been developed. As a result, the OGI detected from various 
approaches differs as much as ±60 days (White et aL 2009). 
Fourth, the trend calculated from the average of phonologic 
metrics over a climate region is not exactly the same as that 
calculated from pixel-based metrics. Although the direction of 
the overall trend in a climate region follows the pixel-based 
trends in majority pixels, the statistical significance and mag- 
nitude of the overall trend are contributed from all pixels. As a 
result, the significance of overall trend in a climate region may 
not match well with significant pixel-based trends. With this 
in mind, we could understand better the discrepancy of phe- 
nological trends and interannual variations from various 
studies. 

This study revealed that the early trend and magnitude of 
rate in OGI varied greatly with climate regions in middle-high 
latitudes. In polar climate, OGI came early and GSL was 
prolonged with significant trends in NA from 1982 to 2010, 
but no significant trends were found in Eurasia. In cold and 
temperate climates, there were more pixels with significantly 
early trends than those with late trends from 1 982 to 1 999, and 


early trends reduced from 2001 to 2010. It is likely that the 
impacts of warming climate on greenness increase were more 
consistent in polar climate than other climate regions. 
Moreover, our results in northern middle-high latitudes are 
generally comparable with some previous studies, but not with 
others. Based on AVHRR GIMMS from 1982 to 1999, the 
early shift rate in OGI was 0.36 days/year in Eurasia and 
0.42 days/year in NA (Zhou et al. 2001), 0.25 days/year from 
45° to 75°N (Tucker et al. 2001), 0.74 days/year in North 
China (Piao et al. 2006), and 0.31 days/year in the temperate 
Northern Hemisphere (de Jong et al. 2011). According to 
AVHRR PAL data, the early shift rate in OGI was 
0.62 days/year in NA and 0.42 days/year in Europe from 
1985 to 1999 (de Beurs and Henebry 2005a), 0.03 days/year 
in temperate climate from 2000 to 2008 (de Jong et al. 2011), 
and 0.54 days/year in Europe from 1982 to 2000 (Stockii and 
Vidale 2004). Based on AVHHR GVIX (4 km) data, in NA 
from 1982 to 2005, OGI came early with a rate of 0.32 days/ 
year from 40°N northward but delayed 0.24 days/year in 
31.5°N southward (Zhang et al. 2007). In contrast, significant 
OGI trends were generally not found across NA according to 
biweekly AVHRR GIMMS data (8 km) from 1982 to 2003 
(Reed 2006) and from 1982 to 2006 (White et al. 2009). 
Similarly, GSL derived from AVHRR GIMMS data showed 
an increased rate of 0.95 days/year in Eurasia and 0.63 days/ 
year in NA (Zhou et al. 2001), 0.05 days/year from 45° to 
75°N from 1982 to 1999 (Tucker et al. 2001), and 0.96 days/ 
year in Europe from 1 982 to 2000 (Stockii and Vidale 2004). 

In the dry climate of western NA, OGI, overall, showed no 
significant trend but GSL increased during 1982-1999. This is 
likely associated with the complex cycles in precipitation and 
dry season length during past decades (Groisman and Knight 
2008; Zhang et al. 20 1 0). In dry climate across Asia, GSL was 
significantly prolonged in BS and BW during the period of 
2000-2010 although significant shortening occurred in BS 
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from 1 982 to 1999. The positive trend is likely associated with 
the precipitation increase during the recent decade on the 
Tibetan Plateau (Piao et al. 2006; Shen et al. 2011). 
Moreover, our results also show that OGI presented no sig- 
nificant trends from 2000 to 201 0 and that OGI shifted early in 
BW from 1982 to 1999. This complex pattern is comparable 
with previous results from GIMMS data which indicates that 
an early GIMMS OGI trend (0.38 days/year) appeared from 
1982 to 1999 (ShenetaL 2011). 

In South America, OGI showed a late trend during the last 
3 decades particularly from 1982 to 1999. Correspondingly, 
GSL decreased 1-15 days in tropical and temperate climates 
from the 1990s to the 2000s although it increased from the 
1980s to the 1990s. This reduction during the 2000s is con- 
sistent with the finding of NPP (net primary productivity) 
reduction because of drought impacts induced from warming 
temperature (Zhao and Running 2010). 

The early OGI trend and prolonged GSL in the Northern 
Hemisphere of Africa were evident from 1982 to 1999; how- 
ever, the rate reduced or converted to the opposite directions 
during 2000-2010. This pattern is generally consistent with 
previous studies (Olsson et al. 2005; Heumann et al. 2007; 
Vrieling et al. 2013). The trend is likely related to the increase 
of rainfall during 1980s to the 1990s as well as land use 
change (Olsson et al. 2005). But the result is discrepant with 
that from GIMMS AVHRR NDVI during 1982-2003. The 
latter indicates that shortened GSL was dominant from 1982 
to 2003 (Julien and Sobrino 2009). 

The overall insignificant GSL change (but more pixels with 
significant increases than decreases) in the Southern 
Hemisphere of Africa from 2000 to 2010 is discrepant with 
the decreased trend of NPP as described by Zhao and Running 
(20 1 0). However, the insignificant trend from 1 982 to 20 1 0 in 
this study is similar to the results of GIMSS data during 1 982- 
2011 (Vrieling et al. 2013) and 1982-2006 (de Jong et al. 
2011) where pixels with significantly increased and decreased 
trends were sparsely distributed across the region. 


Conclusions 

We investigated interannual variations and trends of LSP from 
1982 to 2010 over individual climate regions globally. The 
seasonal vegetative greenness was simulated from long-term 
EVI2 time series after snow and cloud effects were explicitly 
reduced. Phenological metrics were retrieved from PLM- 
LSPD which provides a flexible, repeatable, and realistic 
means to monitor seasonal and interannual dynamics in veg- 
etation across the globe. Because frequent cloud cover during 
a vegetation growing season strongly limits the reconstruction 
of a vegetation growing cycle, the proportion of good obser- 
vations provides the overall reliability of phenology detection 
in individual pixels. On the other hand, the quality of 


phenology retrieval could vary with ecosystem types. Since 
the annual EVI2 cycle is more complex and less distinguish- 
able in drylands and evergreen forests than that in other 
ecosystems, the retrieved phenological metrics likely contain 
relatively high uncertainty. Thus, it is urgently needed to 
develop robust methods to validate and assess the accuracy 
of long-term LSP detections and the phenological trends 
across various ecosystems. 

The spatial patterns in OGI and GSL were very similar 
during the last 3 decades because they were mainly controlled 
by global climate patterns. However, interannual variation in 
phenology was larger in tropical and dry climates than in cold 
and polar climate regions. This is likely attributed to the 
difference in the seasonal variation between temperature and 
precipitation. Temperature seasonality is generally regular 
interannually, which leads to the low variability in vegetation 
phenology in middle-high latitudes. In contrast, rainfall sea- 
sonality is relatively more variable, which results in the high 
variability of interannual OGI and GSL in dry climate and 
seasonal dry tropical regions. 

The phenological trend and interannual variation exhibited 
great variability across geographic and climatic regions. OGI 
generally shifted early and GSL was prolonged from 1982 to 
2010 in temperate, cold, and polar climates in the North 
Hemisphere although the trends were not necessarily signifi- 
cant in individual climate regions and the proportion of pixels 
with significant trends was generally less than 10 %. 
However, the pixels with significantly advanced OGI reduced 
during 2000-2010 relative to those from 1982 to 1999, and 
the significant trends were more evident in North America and 
Asia than those in Europe. In precipitation-controlled dry 
climate over the Northern Hemisphere, OGI and GSL pre- 
sented no significant trends in Europe and western NA (except 
increased GSL in BW of NA) while GSL was mostly 
prolonged in Asia from 2000 to 2010. In South America, late 
OGI was consistent, particularly, from 1982 to 1999 while 
GSL varied between the periods of 1982-1999 and 2000- 
2010, where GSL significantly reduced in temperate climate 
but increased in dry climate from 2000 to 2010. No significant 
OGI trend was found, but prolonged GSL was evident over 
individual climate regions in the Northern Hemisphere of 
Africa during the last 3 decades because there were more 
pixels with significantly early OGI and increased GSL than 
those with late OGI and shortened GSL. OGI in the Southern 
Hemisphere of Africa mainly showed late trends, and the 
proportion of pixels with reduced GSL trends was larger than 
that for prolonged trends from 1982 to 1999, but it was 
reversed from 2000 to 2010. In Australia, GSL exhibited 
considerable interannual variation, but the consistent trend 
lacked presence in most regions. 

The combination of the overall trend and pixel-based trend 
provides phenological characteristics varying at different spa- 
tial scales in a climate region. This could improve our 
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understanding of the phenological response to climate change. 
However, future research is needed to explore how trends 
(including significance, magnitude, and direction) in individ- 
ual pixels affect the overall trend in a climate region and why 
positive and negative trends coexist in the same climate 
region. 

Moreover, the trend was conservatively analyzed by sepa- 
rating the two periods of 1982-1999 and 2000-2010 because 
of the different qualities and some inconsistencies between 
AVHRR LTDR and MODIS CMG. In this case, the direction 
of the overall trend in a climate region from 1982 to 2010 is 
suggested to follow the direction that consistently occurred in 
both periods. However, it was uncertain if the trend direction 
was opposite during these two periods. As a result, farther 
evaluation and calibration of the EVI2 time series derived 
from AVHRR LTDR and MODIS CMG could enable us to 
provide more reliable long-term phenological variation and 
trends. 

Finally, the trends of vegetation phenology during the last 3 
decades were substantially variable and inconsistent in most 
parts of the globe. To understand the complex responses of 
long-term phenology to climate change, it is our ongoing 
efforts to examine the correlation between phenologic metrics 
and climate variables at individual climate regions globally. 
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